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I.  THE  PARABOLIC  WAVE  EQUATION 


The  wave  equation  for  the  complex  acoustic  pressure  field  p  away  from 
acoustic  sources  is 


(1) 


where  C  is  the  speed  of  sound. 

2  2  2 

Here  p  =  p(x,y,z,t) ,  C  =  C(x,y,z) ,  and  V2  =  —  +  —  +  —  . 

3xz  3yZ  3zZ 

We  use  the  coordinate  system  shown  below,  in  which  z  is  the  horizontal  (range) 
coordinate  and  x  is  the  vertical  coordinate: 


Since  p  is  assumed  complex,  the  observed  signal  is  Re(p),  the  actual  pressure. 
Following  Dashen  [1]  and  others,  we  assume  that 

p(x,y,z,t)  =  ^(z,x,y)el(-kz"a)t'1  ,  (2) 

where  ip  is  a  complex  envelope,  0)  is  the  frequency,  and  k  =  1/A  is  the  wave 
number,  A  being  wavelength. 

Letting  Aeo  =  oo/k  =  C^,  a  reference  sound  speed,  we  define  the  index  of 
refraction  n  by 

n(x,y,z)  =  CQ/C(x,y,z)  . 

Then  (1)  becomes 

V2p  +  k2n2p  =0,  (3) 


the  reduced  wave  equation. 
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2  2 
d  a 

Putting  A  =  — ^  j  and  using  (2) ,  we  see  that  (3)  may  be  written 

3xZ  9yZ 


+  -^-1  +  2ik  ||  +  k2  (n2-'l)  \p  =  0 


Aijj  •  2 

9z 


Suppose  we  make  the  following  assumptions: 


2  2  9  f 

A)  The  term  9  <Jj/9z  may  be  neglected  in  comparison  with  2ik  ,  i.< 


d2lp 

9 

« 

2ikS 

9z2 

and 


B)  ^(n2-l)  s  n-1  ,  i.e.,  n  =  1 . 


Thus,  we  have 

1  If +  A +  kcn-i>'j’=o > 

or,  defining  y  by 

b(x,y,z)  =  l-n(x,y,z),  (4) 

we  get 


Equation  (5)  is  called  the  parabolic  wave  equation.  As  we  have  seen,  approxi¬ 
mating  the  reduced  wave  equation  (3)  by  the  parabolic  wave  equation  (5)  rests 
on  the  assumptions  A)  and  B)  above.  Assumption  B)  says  just  that  the  sound 
speed  varies  little  throughout  the  medium,  i.e.  the  approximation  C  = 
introduces  a  small  percentage  error  everywhere  in  the  medium.  This  assumption 
is  surely  valid  in  underwater  acoustics.  For  example,  if  one  considers  the 
representative  speed  of  sound  profile  on  p.400  of  [2],  one  easily  estimates 
that  the  approximation 


t 
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C  a  C  . 
mm 


=  1475  m/sec. 


introduces  an  error  of  less  than  5%  down  to  a  depth  of  6  km. 

The  validity  of  the  parabolic  approximation,  and  of  assumption  A)  in 
particular,  has  been  carefully  considered  by  Tappert  in  the  survey  article 
[3],  and  by  Flatte  et  al  [4].  In  the  latter  reference,  pp.  76-77,  it  is 
stated  that,  for  typical  conditions  for  sound  propagation  in  the  ocean. 
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so  that  assumption  A)  seems  quite  safe  in  underwater  acoustics.  Tappert 

observes  that  the  parabolic  approximation  applies  essentially  to  long-range 

(R  =  range  »  A,  typically  R  >  50-60  km.)  propagation  in  the  ocean  sound 

channel,  which  acts  like  a  vertically  thin,  horizontally  elongated  waveguide; 

consequently  if  we  think  in  terms  of  rays,  the  angle  9  between  the  ray 

direction  and  the  horizontal  is  small  (0  <  16°)  for  rays  which  are  not 

max 

attenuated  to  zero  intensity  over  the  long-range  propagation.  The  envelope 
ip  varies  more  slowly  with  the  range  coordinate  z  than  with  depth  x.  The 
parabolic  approximation  retains  diffraction  effects,  unlike  geometric  acoustics, 
which  assumes  that  wavelengths  are  small  enough  that  diffraction  can  be  ignored 
(see  below).  The  parabolic  approximation  is  therefore  valid  to  significantly 
lower  frequencies  than  geometric  acoustics.  Finally,  the  parabolic  approxi¬ 
mation  is  not  limited  to  stratified  media,  as  are  normal  mode  expansions. 

II.  GEOMETRIC  ACOUSTICS  APPROXIMATION  TO  THE  PARABOLIC  WAVE  EQUATION 

Since  the  parabolic  wave  equation  retains  the  diffraction  effects  whose 
neglect  is  the  essence  of  the  geometric  acoustics  approximation  (at  least 
in  the  presence  of  gentle  refraction,  n  s  as  we  have  assumed  -  see  [2], 
P-118),  we  should  be  able  to  derive  geometric  acoustics  from  the  parabolic 
wave  equation  (5)  via  a  suitable  neglect  of  diffractive  terms.  We  proceed 
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to  do  so,  following  Tappert  [3,  p.234  ff.]. 
Let  us  write 


^(z,x,y)  =  A(z,x,y)e10(-Z,X,y')  ,  (6) 

where  A  is  a  real  amplitude  and  0  is  a  real  phase.  For  simplicity,  here  and  in 
the  sequel  we  ignore  the  (azimuthal)  dependence  of  iJ>,A,<j>,  etc.  on  y:  we  assume 

ip  =  \p(z,xj  =  A(z ,  A,  0  real  .  (7) 


Now  substitute  (7)  into  (5),  and  divide  by  e^j  the  real  part.bf  the  resulting 
equation  yields 


-  ^ 


pkA  =  0 


(8) 


and  the  imaginary  part  yields 

3A  .}_  3A  _3£  A  3 V  „ 

3z  +  k  9x  3x  +  2k  „  2 
3x 


(9) 


Now  define  0  by 

0(z,x)  =  1 .  (10) 

0. 

0  .  (H) 

Also,  in  terms  of  0,  (8)  is 

,  S*  1  32A  1,  ..2  . 

-A  it  *  a  77  -  5kA0  - vkA  ■ 0  • 

3x 

If  we  differentiate  this  with  respect  to  x,  divide  by  (-kA) ,  and  use  (8) 
again,  we  obtain 

30  19  -  3d  J_  3. 

3z  +  3x  3x  +  2  3x 

2k 


I 

A  3x2 


(12) 


Then  (9)  becomes 

3A  fl  3A  A  30 

3z  +  3x  +  2  3x 

Multiplying  this  by  2A,  we  obtain 


3  ,..2.  3  ,a.2. 

Si  (A  )  *  3-  (8A  3  - 
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Now  the  geometric  acoustics  approximation  is  that  k  »  1,  and  so  the  second, 
diffractive  term  on  the  right-hand  side  of  (12)  may  be  dropped: 

§ +  6  H =  -  S  •  ^i2GE^ 


Consider  the  curves  x  =  x(z)  defined  by 

dx 


dz 


0(z,x)  . 


(13) 


Along  these  curves  (12GE)  says 

d 


or  by  (13)  again. 


dz 


0(z,x(z))  =  -  3- 


9p 

9x  » 


d  X  _  -9y 


dz 


2  9x  * 


(14) 


which  is  the  differential  equation  for  the  rays  of  geometrical  acoustics  in  a 
horizontally  stratified  ocean  (i.e.,  n=n(x))  for  which  n  =  1  and  0  <  1,  as 
we  assume  for  the  parabolic  approximation.  Before  verifying  this  claim,  we 
notice  that  the  geometrical  acoustics  assumption  k  »  1  in  (10)  implies  that 
0  <  1,  so  that  tan  0=0;  hence  (13)  says  that  0  is  just  the  angle  between 
a  ray  and  the  horizontal. 

In  general,  the  ray  equations  are  (see,  e.g.,  [5],  p.79) 


d  ,  dx. 

n  Si  (n  15> 


,  3 

2  9x 


(n  ), 


(15) 


2  2  2 

assuming  n  =  n(x),  where  a  is  arclength  along  a  ray.  Since  (da)  *(dz)  +(dx)  , 


and  4^  =  tan  0,  we  have 
dz 


r  j  >  2 

da 


dz 


1  + 


dx 


dz 


2  2 

=  1  +  tan  0  =  sec  0  , 


so 


dz 

da 


cos  0,  and  therefore 


dx  _  dx  dz_  dx 
da  ~  "dz  da  “  dz 


cos  0 
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Thus,  (15)  becomes 


n 


da 


a  dx 
n  cos  9  t— 
dz 


i  3  ,  2, 

" 2  s (n  ’ 


But  Snells  law  is 


s  =  n  cos  0  =  constant , 
so  the  above  is  simply 


ns 


d 

do 


dx 

dz 


> 


or 


cEx 

dz2 


_  1 
-  2 


3  f  2, 


2 

But  n  =  1  and  0  1  imply  s  =  1,  and  l(n  -1)  =  n-1  =  -y,  and  so  the  above 

equation  reduces  to  (14),  as  claimed  above. 

Before  leaving  the  geometric  acoustics  approximation  temporarily,  notice 

2 

that  (11)  is  just  a  transport  equation  for  the  acoustic  power  A  along  the 
rays . 


III.  RANDOM  RAY  MECHANICS 

We  are  going  to  propose  a  theory  of  rays  z  x(z)  in  which  x(z)  is  a 
random  process  indexed  by  the  range  coordinate  z.  This  theory  is  Nelson* s 
stochastic  mechanics  [e.g.,  6-9],  and  the  presentation  initially  follows 
Nelson  [8] .  Our  goal  is  to  recapture  diffractive  effects  for  geometric 
acoustics  by  regarding  x(z),  for  fixed  z,  as  spread  out  vertically  according 
to  a  probability  density  p(x,z).  We  shall  see  mathematically  that  such  a 
picture  is  equivalent  to  the  parabolic  wave  equation  (5).  Of  course.  Nelson's 
theory  is  a  stochastic  version  of  the  Schrodinger  equation  of  quantum 
mechanics,  but  this  equation  is  formally  identical  to  (5),  provided  we  replace 
the  time  coordinate  by  the  range  coordinate  z.  We  shall  omit  most  proofs  in 
this  section,  referring  the  interested  reader  to  [6-9] . 
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Consider  a  random  path,  or  ray,  governed  by  the  Ito  stochastic  differential 
equation  [10-12] 


dx(z)  =  b(x,z)dz  +  /2v  dw(z). 


(16) 


where  b  is  called  the  drift,  v  is  a  constant  called  the  diffusion  coefficient 
(we  retain  these  terms  from  the  standard  literature  where  z  is  replaced  by  t, 
although  calling  b  the  "spread"  and  v  the  "diffraction  coefficient"  would 
perhaps  be  more  appropriate  here),  and  w(z)  is  standard  Brownian  motion. 

Let  p(x,z)  be  the  probability  density  for  x(z) .  Then  p  satisfies  the  equation 


(17) 


which  physicists  call  the  Fokker-Planck  equation  and  mathematicians  call 
Kolmogorov's  forward  equation  (see,  e.g.  [13],  p.275). 

We  may  also  describe  the  process  x(z)  with  the  direction  of  z  reversed 
(analogous  to  time-reversal  when  z  is  replaced  by  t) ;  since  the  reversed 
process  is  again  of  the  same  type  (cf.  [13],  p.83),  we  have 


dx  =  b*(x,z)dz  +  /2v  dw*  , 


(18) 


for  a  suitable  backward  drift  b*  and  a  (possibly  different)  Brownian  motion 
w*,  where  we  understand  dx  =  x(z)  -  x(z-dz),  dz  >  0.  We  then  have  the  backward 
Fokker-Planck  equation 


(19) 


Let  v,  the  current  velocity,  be  defined  by 


v  =  (b+bJ/2  . 


(20) 


Then  (17)  and  (19)  yield  the  equation  of  continuity, 


(21) 


Let  the  osmotic  velocity  u  be  defined  by 


u  =  (b-bJ/2 


(22) 
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Then  [7,  p.105]  we  have 


v  _3p_ 

“  p  3x  ’ 


(23) 


The  formula  (23)  for  Brownian  motion  was  in  essence  derived  by  Einstein 

([14],  p.9  ff.);  the  term  "osmotic  velocity"  is  used  because,  in  the  case  of 

Brownian  motion  of  particles  suspended  in  a  liquid,  u  is  the  velocity  particles 
must  acquire  to  overcome  osmotic  effects  ([7],  p.21). 


We  now  assume  that  there  exists  a  function  V(x,z)  such  that 


^  Pti^+lv2  +  V)p(x,z)dx  =  0  •  (24) 

-00 

In  the  diffusion  case  (z  -*■  t)  ,  V  is  the  potential  energy,  and  (24)  expresses 
energy  conservation:  the  average  total  energy  is  constant  in  time.  We  shall 
see  that  in  our  case  we  may  take  V  =  y,  the  variation  y  =  1-n  of  the  index 


of  refraction  from  its  typical  value  1  (recall  that  n  a  1  in  the  parabolic 
approximation).  From  (24)  it  follows  ([8],  p.174)  that 


,  32u  3u.  ,3v  3v,  3V  n 

-Cv  —  +  u  ■§£■)  +  (—  +  v  — )  +  —  =  0 


3x 


'3z 


3x 


3x 


(25) 


If  we  differentiate  (23)  with  respect  to  z  and  then  use  (21) ,  we  find  that 


3u 

3z 


3  v  3  ,  v 

*  -« —i  -  s (uv)  ■ 


3x 


Assuming  that  p  never  vanishes,  let  R  be  given  by 

R  =  2  log  P, 

so  that  (23)  becomes 


u  =  2  v 


3R 
3x  ’ 


Also  write  v  as  a  gradient:  choose  a  functions  such  that 

3  S 


v  =  2  v 


3x  * 


(26) 

(27) 

(28) 

(29) 


(Note:  the  analogues  of  (28)  and  (29)  in  [8]  are  unfortunately  in  error.) 
Now  define  ip  by 


\p  =  e 


R+iS 


(30) 
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and  choose 


(31) 


Then  a  simple  calculation  shows  that  ip  satisfies  the  parabolic  wave  equation 
(5): 


,  ,  1  djjj 

9z  2k  .  2 

dX 

provided  we  take  V  to  be 


-  ykip  =  0  , 


(5) 


V  =  y 


(32) 


(Actually,  since  we  must  integrate  (25)  and  (26)  with  respect  to  x  to  derive 
(5),  a  seemingly  arbitrary  complex  constant  of  integration  a  =  a(z)  is 
introduced,  so  that  we  get 

(*)  i||  +  -Mr  -  Ukip  +  a(z)  \p  =  0 

3x 

—  000 

But,  because  p  =  #  is  a  probability  density,  J  p(z,x)dx  =  1,  for  all  z. 

-00 

Since  (*)  multiplied  by  ip  and  its  complex  conjuate  multiplied  by  ip  imply 


.  9p  _  1  3i 

1  3z  "  “  2K  3x 


+  (a-a)p  , 


we  have 


C)  *oo  —  — — 

0  =  i9- J  p  dx  =  (a-a) J  p  dx  =  (a-a)  , 


_  #oo 


so  a  must  be  real.  But  S  in  (29)  is  determined  only  up  to  an  arbitrary  additive 
function  of  z.  Choosing  S  appropriately,  we  may  insure  a  =  0.  This  argument  is 
due  to  Nelson  [7,  p.131].) 

Conversely,  suppose  ip  satisfies  (5),  and  define  R,  S,  u,  v,  by  (28),  (29), 
and  (30).  Then  u  and  v  satisfy  (25)  and  (26),  with  V  =  y.  Thus,  each  solution 
x(t)  of  (16)  gives  rise  to  a  solution  \p  of  the  parabolic  wave  equation  (5), 
and  conversely:  knowing  ip,  we  know  u  and  v,  and  hence  b.  In  this  sense,  (5) 
is  strictly  equivalent  to  a  theory  of  random  rays  governed  by  (16) .  From 
the  point  of  view  of  random  ray  theory,  solving  (5)  for  ip  is  merely  a 
convenient  method  of  determining  the  drift  b,  and  hence  what  stochastic 
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differential  equation  the  rays  obey.  Solving  the  linear  equation  (5)  is 
simpler  than  solving  directly  the  coupled,  non-linear  system  (25),  (26)  of 
partial  differential  equations. 

It  will  be  convenient  to  have  an  alternative  account  of  the  dynamics 
embodied  in  (25).  IF  f(x,z)  is  a  smooth  function,  define  the  mean  forward 
derivative  Df(x(z),z)  of  f(x(z),z)  by 


Df(x(z),z)  =  lim  E  lf(x(z+Az),z+Az)-f(x(z),z) 
Az+0 


Az 


(33) 


and  the  mean  backward  derivative  D*f(x(z),z)  of  f(x(z),z)  by 


D„f(x(z),z)  =  lim  (  ftxW-z)  -  f(x(z-Az),z-tz) 
AzlO  ^  Az 


Fz  ' 


(34) 


where  P  =  a{  x(s)  1 0<s<z}  is  the  a-field  generated  by  x(s)  for  s  <.  z,  F  = 
a{x(0)|s>z}  is  the  a-field  generated  by  x(s)  for  s  >.  z,  and  E{X | F }  is  the 
conditional  expectation  of  the  random  variable  X  with  respect  to  the  a-field  F. 
Now,  by  Ito's  lemma  (e.g.,  [11],  p.33),  we  have 

2 

Df(x(z),z)  =-||(x(z),z)  +  b(x(z),z)  ||(x(z),z)  +  V  — j(x  (z)  ,  z)  ,  (55) 

9x 

and  similarly, 

2 

D*f(x(z),z)=  ||(x(z),z)  +  b*(x(z),z)  ||(x(z),z)  -  v-^-|(x(z)  ,z)  .  (36) 

9x 

In  particular,  if  f(x(z),z)  =  x(z),  then 

Dx(z)  =  b(x(z) ,z)  (37) 

and 

D*x(z)  =  b*(x(z),z)  .  (38) 

Now,  define  the  stochastic  acceleration  a  by 


a  =  i{DD*+D*D}(x(t)) • 


(39) 
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Then  by  (37)  and  (38) , 


a  =  J{Db*  +  D*b}  , 


or,  by  (35)  and  (36), 


2  2 

In  terms  of  u.  =  (b-b*)/2  and  v  =  (b+b*)/2,  bb*  =  v  -u  ,  so  this  is 

? 


(40) 


Hence,  (25)  is  simply 


a  3x  ’ 

or,  defining  the  force  F  by  F  =  -9V/9x  , 


(41) 


F  =  a 


(42) 


which  is  just  Newton's  second  law  for  a  particle  of  unit  mass.  Thus,  with 
the  stochastic  acceleration  defined  by  (39) ,  the  dynamical  equation  (25)  is 
simply  Newton's  second  law,  rather  than  energy  conservation  (eq.  (24)). 

We  also  have  an  important  formula  relating  the  ordinary  z-derivative 
of  average  quantities  to  averages  of  their  mean  forward  and  mean  backward 
derivatives ( [7] ,  p.98;  [9],  p.204):  if  f(x,z)  and  g(x,z)  are  smooth  func¬ 
tions,  and  E  denotes  averaging  in  the  probability  space  of  the  process  x(z), 


then 


In  particular,  with  g  =  1  this  yields 


^  E(f(x(z) ,z))  =  E(Df(x(z) ,z))  , 


(44) 


and,  interchanging  f  and  g  and  putting  g  =  1,  we  find 


^E(f(x(z),z))  =  E(D*f(x(z),z)). 


(45) 


Simple  consequences  of  (44)  and  (45)  are 
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E(f (x(z) , z))  =  E(l(D+D*)£(x(zD,z)) 


(46) 


and 


E((D-Djf(x(z),z))  =  0. 


(47) 


In  particular,  for  f(x(z),z)  =  x(z),  (46)  and  (47)  imply,  by  (37)  and  (38), 


E(x(z))  =  E(v(x(z) ,z)) 


(48) 


and 


E(u(x(z)  ,z))  =  0. 


(49) 


Differentiating  (48)  with  respect  to  z,  and  using  (46),  we  obtain 

.2 


dz 


E(x(z))  =  E(i(D+D*)v(xCz),z))  , 


or,  by  (35)  and  (36), 


~2  E(x(z))  =  E(|j(x(z),z)  +  v(x(z),z)  |^(x(z),z)).  (50) 

dz 

Now,  let  us  take  expectations  in  (40)  to  get 

E(a(z))  =  E(|j(x(z),z)  +  v(x(z),z)  |^(x(z),z))-  E(u(x(z)  ,z)|£(x(z)  ,z) 

~2 

+  v  •  (x (z)  ,  z) )  . 

8x 

We  shall  show  that  the  average  in  the  second  term  is  zero.  Indeed,  by  (47) 
with  f  =  u, 

0  =  E(i(D-D*)u(x(z),z)), 


and  by  (35)  and  (36) ,  this  is 

2 

0  =  E(u(x(z) ,  z)  |^(x(z),z)  +  v  —•  (x(z),z)), 

Bx 

as  claimed.  Hence 

E(a(z))  =  E(|^(x(z),z)  +  v(x(z) ,z)|^(x(z) ,z)) 

Comparing  (50)  and  (51) ,  we  have 

,2 


dz" 


E(x(z))  =  E(a(z))  , 


(51) 
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or,  recalling  (41), 

,2 


dz 


2  E(x(z))  = 


-E(|£(x(z),z)) 


(52) 


This  equation  says  that  the  average  ray  z  ->  E(x(z))  feels  the  average  force, 
-  9V 

and  not  -^-(E(x(z))  ,z) ,  the  force  along  the  average  ray. 


IV.  RANDOM  RAYS  AND  THE  GEOMETRIC  ACOUSTICS  APPROXIMATION 

Suppose  we  have  a  solution  of  the  parabolic  wave  equation  (5).  Thus, 
we  have  vector  fields  u  and  v,  defined  by  (28) -(30),  so  we  may  define  three 
different  kinds  of  paths,  or  rays: 


G  G 

geometric  acoustics  rays  x  :  z  x  (z)  : 

These  are  defined  as  solutions  to  the  ray  equation  (14) 

d2xG,  ,  -3p,  G 

— T(z)  =  ^(x  (z)). 

dz 


(53) 


P  P 

parabolic  rays  x  :  z  ■*  x  (z) : 


These  are  defined  by  (13) ,  so  that 

dxP  ,  P,  .  . 

=  v (x  (z) ,  z) 


dz 

(note  that  0  of  section  II  is  the  same  as  v  of  section  III:  cf.  (6)  and 
(10)  with  (29) -(31)),  and  thus,  by  (12),  satisfy 


(54) 


d  x  ,  .  3p,  P,  ,,  13 

— Cz)  =  -  ^(x  (z))  + -2  R 


dz 


2k 


1  3  2 A1 


A  .  2 

dX 


(x  (z) ,z)  , 


or, 

d2  r  .  3y,  P,  .  1  9^u,  P,  .  .  ,  P,  ,  .3u,  P,  ,  .  ,Pr. 

- 2^z)  =  "  2k  — 2^X  (z)’z)  +  u(x  (z)>z)-jw(x  (z)>z)  »  (55) 

dz  3x 

1  9R.  1  9 

where  the  last  equation  follows  from  u  =  ^  (log  A),  by  direct  compu¬ 


tation,  or  by  comparing  (12)  and  (25) . 
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C.  random  rays  x:  z  x(z)  :  defined  as  a  random  process  as  in  section  III, 

in  particular  the  average  random  ray  E(x(z))  satisfies  (52): 

2 

^ECxCz))  =  -E(g-(x(z),z))  .  (56) 

dz 

All  three  objects  xG,  xP,  and  E(x(z))  are  in  general  distinct,  since  they 

satisfy  the  different  equations  (53) ,  (55) ,  and  (56) . 

Since  the  parabolic  wave  equation  is  just  the  Schrodinger  equation  of 

quantum  mechanics,  and  since  we  expect  ray  theory  to  be  meaningful  when  k  >>1, 

i.e.,  when  1/k  «  1,  where  (1/k)  is  playing  the  role  of  D  =  Planck's  constant 

divided  by  2ir  in  (5) ,  we  can  interpret  some  approximate  relationships  among 

XP,  x(z),  and  xG  in  terms  of  approximating  quantum  mechanics  by  classical 

mechanics  when  "fi  «  1.  For  example,  the  approximation 

P  ..  G 
x  =  x 

achieved  as  in  section  II  by  dropping  the  last  two  terms  on  the  right  side 
of  (55)  has  been  well  studied  in  the  quantum-mechanical  context.  It  can  be 
understood  as  ignoring  a  "quantum-mechanical  potential,"  thereby  converting 
a  non-linear  fluid-mechanical  system  (the  "Madelung  fluid")  into  a  linear 
fluid-mechanical  system  (the  "Hamilton-Jacobi  fluid")  whose  dynamics  consti¬ 
tute  a  well-known  version  of  classical  mechanics  (see,  e.g.,  [15]  and 
references  therein) . 

2 

G  Id 

The  approximation  x(z)  =  x  (z)  amounts  to  replacing  y(DD*  +D*D)  by  —  ^  , 

i.e.,  replacing  stochastic  acceleration  by  the  usual  second  time  derivative. 
Studying  the  semiclassical  limit  of  quantum  mechanics  in  its  stochastic 
mechanical  version  is  a  relatively  untouched  field  of  research  (however,  see 
the  interesting  paper  [16]). 

Q 

One  might  hope  that  E  x(z)  =  x  for  appropriate  initial  conditions,  but 
this  is  false  in  general:  the  average  random  ray  is  not  a  geometric  acoustics 
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ray,  since  they  satisfy  different  differential  equations,  (53)  and  (56) .  The 
failure  of  E(x(z))  to  be  a  geometric  acoustics  ray  is  measured  by 

E(|H  (x(z),z))  -  |H(B(x(z),z))  , 

the  error  made  in  approximating  the  mean  force  by  the  force  on  the  mean  ray. 

Estimation  of  the  size  of  this  error  seems  worth  further  investigation,  as  do 

P  r  G 

the  approximation  x  s  x  and  x(z)  =  x  ,  understood  as  semiclassical 
quantum-mechanical  limits. 

A  related  question  is  the  estimation  of 

E(x(z))  -  xP(z), 

P 

which  might  well  proceed  from  the  formula  (recall  that  x  (z)  is  non-random) 

Z  P 

E(x(z)-xP(z))  =  /  E[v(x(s),s)-v(xP(s),s)]ds  , 
s=0 

together  with  the  assumption  that  v  is  Lipschitz. 
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V.  THE  FEYNMAN- KAC  FORMULA  AND  SOLUTIONS  TO  THE  PARABOLIC  WAVE  EQUATION 

Invented  in  1948,  Feynman's  path  integral  [17,  18]  has  long  since 
established  itself  as  a  powerful  alternative  formulation  of  quantum  dynamics, 
embodying  direct  heuristic  appeal  and  achieving  success  in  an  impressively 
broad  range  of  applications  [see,  e.g.,  19].  Indeed,  Dashen's  pioneering 
application  of  Feynman's  path  integral  to  wave  propagation  in  random  media  [1] 
is  the  primary  stimulus  of  the  present  work.  As  Dashen  has  shown,  the 
Feynman  integral  representation  of  the  solution  of  the  Parabolic  Wave 
Equation  with  random  fluctuation  y  allows  one  to  calculate  moments,  corre¬ 
lation  functions,  and  various  other  important  characteristics  of  the  random 
field  ty. 

From  a  mathematical  point  of  view,  however,  Feynman's  path  integral  has 
proved  profoundly  problematical.  It  was  soon  realized  [20]  that  Feynman's 
integral  cannot  be  understood  as  an  ordinary  integral,  since  it  arises  from 
no  measure.  Various  attempts  at  rigorous  formulations  of  Feynman  integrals 
via  generalizations  of  the  concept  of  measure  have  proved  fertile  fields  for 
mathematical  research.  But  the  most  significant  development  for  our  purposes 
is  Kac's  discovery  [21;  22,  pp.  165-172]  that,  provided  we  begin  with  the  non- 
homogenous  heat  equation  instead  of  the  parabolic  wave  equation,  Feynman's 
procedure  does  lead  to  a  well-defined  integral,  the  Wiener  integral,  arising 
from  a  well-defined  measure,  Wiener  measure.  The  resulting  formula  for  the 
solution  of  the  heat  equation  is  called  the  Feynman-Kac  formula  (equation  (FK) 
below) .  This  history  is  summarized  succinctly  in  the  Introduction  to  Barry 
Simon's  admirable  book  [23],  where  many  additional  references  are  provided. 

Since  the  passage  from  the  parabolic  wave  equation  to  the  heat  equation 
formally  amounts  to  replacing  z  by  -iz  (i.e.,  in  the  quantum- mechanical 
setting,  replacing  t  by  -it),  methods  relying  on  the  Feynman-Kac  formula 
are  often  called  imaginary-time  methods.  A  currently  flourishing  school  of 
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research  in  quantum  field  theory  is  based  on  this  imaginary-time  approach 
([24];  see  also  Nelson's  book  review  [25]).  A  maneuver  frequently  used  by 
those  who  think  of  their  work  as  involving  imaginary  time  is  analytic  continua¬ 
tion  from  an  imaginary  half-axis  to  a  real  half-axis.  We  make  no  use  of 
imaginary  time  per  se  or  of  analytic  continuation  below.  Instead,  we  begin 
with  the  Feynman-Kac  formula  for  the  heat  equation  and,  via  a  trick  of  Donsker 
and  Varadhan  [26],  proceed  along  lines  which,  as  far  as  we  know,  are  new  in 
the  context  of  mathematical  physics. 

Wiener  integral  approaches,  some  of  them  merely  formal  manipulations 
assuming  imaginary  time,  have  been  used  previously  in  the  theory  of  wave 
propagation  in  random  media;  see  the  survey  article  of  Frisch  [27]  for  appli¬ 
cations  up  to  1968.  The  use  of  rigorous  Wiener  integral  techniques,  rather 
than  Feynman  integral  methods,  affords  certain  advantages:  all  of  the  results 
of  the  mathematical  theory  of  measure  and  integration  are  automatically  at 
one's  disposal.  For  example,  as  Chow  [28]  has  pointed  out,  one  may  justify 
the  interchange  of  Wiener  integration  and  expectation  in  the  probability 
space  on  which  the  random  fluctuation  y  is  defined  by  appealing  to  Fubini's 
theorem.  In  any  case,  the  use  of  rigorous  Wiener  integrals  holds  promise 
for  avoiding  the  "balancing  act,  .  .  .  that  is,  our  attempts  to  draw  firm 
conclusions  from  not  so  firm  a  theory"  [19,  p.63],  which  inevitably,  it  seems, 
accompanies  the  use  of  Feynman  integrals.  In  this  spirit,  the  discussion 
that  follows  aims  to  be  rigorous,  modulo  regularity  assumptions  (or  theorems) 
about  the  functions  involved. 

Consider  the  inhomogenous  heat  equation 


3Y  =  1  32Y 
3z  2  3x2 


V(x,z)Y 


(57) 
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where  the  unknown  function  Y  =  Y(x,z)  is  subject  to  the  initial  condition 
that  Y(x,o)  is  a  prescribed  function.  Then  the  Feynman-Kac  formula  for  the 
solution  to  (57)  is 

H  r  '  1  (FK) 


Y(x,z)  =  E^  exp|-/oV(w  (s)  ,s)'dsj-Y(w(z)  ,0) 


where  w(s)  is  Brownian  motion  starting  at  x,  and  "E  "  denotes  integration 

X 

(or  expectation)  in  the  probability  space  on  which  (w(s)}  is  defined,  with 
respect  to  the  corresponding  Wiener  measure.  The  right-hand  side  of  (FK) 
is  thus  a  Wiener  integral  (or  functional  integral”,  or  ”path-space  integral”). 
(We  note  that  V  is  usually  assumed  to  be  independent  of  z;  the  more  general 
case  V  =  V(x,z)  was  considered  by  Faris  [29].) 

Now  suppose  that  u  =  u(x,z)  is  an  ”arbitrary”  function  (modulo  smoothness 
assumptions).  Then,  as  Donsker  and  Varadhan  have  observed  [26,  p . 19 ]  ,  if 
we  put 

j*  — 

r  o 

j 2u(x ,  z)  , 


V (x, z) 


f  u  3u 

9x*  '  92 


(58) 


then  the  solution  to  (57)  is  just  Y(x,z)  =  u(x,z).  Hence,  by  (FK) ,  we  have 


u (x , z)  =  E 


exp{-/0 

f  2  2  1 

9  u/9x  -29u/3z 

(w(s) , s)  dsl  u(w(z)  ,0  ) 

L  1  1 

2u  J 

J 

(59) 


Now,  suppose  also  that  u(x,z)  =  iJj(x,z),  a  solution  to  the  Parabolic  Wave 
Equation  (we  take  k=l  here  for  convenience,  and  continue  to  neglect 
y- dependence) : 

gm  1  8  i I) 

1  Jz  +  2  — 2  "  ^  =  °>  Fi  =  y(x,z).  (60) 

3x 


Then  a  short  computation  shows  that 

g2^ 


3x 


j# 

2  '  2  9z 


/ 


2*  =  -(1*1)^ 


log(ipelM) 


+  ip  , 


(61) 
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where 


so  that 


M(x,z)  =  /Qvi(x,s)  ds, 


3M 


"3z"  =  ^ 


Substituting  (61)  into  (59)  ,  we  find 
^(x,z)  =  E 


(62) 


(63) 


exp{/0 

(l+i)  (log(^eiM))“  iy 

(w(s)  ,s)  dsl  ^(w(z)  ,0) 

) 

In  particular,  for  y  =  0,  the  solution  ^(x,z)  of 

.  3*o  1  3^o 


IT 


2.2 

3x 


=  0 


is  just 


VX’Z)  =  EX 


exp{  /q  [C1+i)-|^-Clog  ^0)]  (w(s)  ,s)ds}iJJo(w(z)  ,0) 


(64) 


(65) 


The  case  of  y  a  Gaussian  random  field,  independent  of  x: 

In  this  case  we  can  use  the  above  approach  to  derive  "the  usual  formula 
obtained  in  the  Markov  approximation"  [1,  p.898]  for  the  first  moment  <\p>  of 
ijj.  Here  and  in  what  follows  "<  >"  denotes  expectation  in  the  probability  space 
on  which  y  is  defined. 


Since  we  assume  here  that  y=y(z) ,  a  short  calculation  shows  that 

^0(x,z)  =  iJj(x,z)elM(-z-) 

(provided  that  ^Q(x,0)  =  iJj(x,0)  ).  Thus  (64)  becomes 


^(x,z)  =  E 


exp{  -i/oP(s)ds}exp{/^[(l+i)^|(log^o)]  (w(s)  ,s) ds}^Q(w(z),0  ) 


or 


iKx,z)  =  exp{-i/oy(s)ds}Ex  exp{/^{(l+i) -^(log^)]  (w(s)  ,s)  ds}^(w(z)  ,0) 


(66) 


(67) 
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because  the  factor  involving  y  is  independent  of  the  path  w (s) .  Using  (65) , 
we  obtain  simply 

ip  (x,z)  =  exp{-i/  y(s) ds}  ip  (x,z).  (68) 

o  ° 


We  assume  that  <y>  =  0  (so  that  by  (4),<f|>  =  1,  the  index  of  refraction  of  a 


homogenous  medium).  Therefore,  /^y(s)ds  is  Gaussian  with  mean  0.  But,  if  Z .is 


a  Gaussian  random  variable  with  zero  mean,  then 

<e-iZ)  - 


-Y  <Z2> 


(69) 


Therefore,  if  we  define  the  function  $  by 


Hz)  =  <(^Zy(s)ds)2>1/2  , 


(70) 


r  ^ 

take  "<>"  of  both  sides  of  (68),  and  use  (69)  for  Z  =  }  y(s)ds,  we  obtain 

o 


(71) 


the  "usual  formula"  mentioned  above  for  the  first  moment  of  the  field  ip  from  the 
theory  of  wave  propagation  in  a  random  medium.  0  is  called  the  rms  phase  fluc¬ 
tuation  (as  computed  in  first-order  geometric  optics)  and  serves  as  a  measure  of 
the  strength  of  the  fluctuation  y. 


The  case  of  more  general  y: 

To  derive  the  formula  (71)  for  more  general  y  =  y(x,z),  it  has  become 
customary  to  make  use  of  a  "Markov  approximation."  According  to  Dashen 
[1,  p.896],  if 

a(((x-x')2  +  (z-z')2)1^2)  =  ^y(x,z)  yi(x*,z'))  (72) 

is  the  correlation  function  of  y,  then  the  Markov  approximation  is  the  replacement 


a (( (x-x* ) 2  + 


a  (|x-x' |) , 


(73) 
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where 

dC|x|)  =  f° a  ((|x|2+z2)1/2)dz  .  (74) 

-00 

In  particular,  for  a  Brownian  path  w(s)  in  this  approximation, 

/y(w(s),s)  y  (w(s'),s')^  =  6(s-s')  8  ( |w(s)-w(s')  j)  =  6(s-s')  6  (0),  (75) 

i.e.,  y's  at  different  points  of  the  same  path  are  uncorrelated  (hence  indepen¬ 
dent,  if  Gaussian),  and  the  above  expression  is  independent  of  w(s)  andw(s'). 
Thus,  for  the  path  integral  representation  (64)  of  y,  equation  (75)  should  amount 
to  ignoring  the  x-dependence  of  y,  so  that  effectively  y=y(z)  as  above.  We 
shall  examine  this  idea,  and  the  Markov  approximation  in  general,  in  detail  in 
the  sequel. 

Important  Paths : 

Elsewhere  [1,  p.898],  Dashen  introduces  the  Markov  approximation  as  follows: 
"The  parabolic  wave  equation  assumes  that  the  normals  to  the  wave  fronts  point 
in  directions  that  are  close  to  the  z-axis.  In  terms  of  the  path  integral  this 
means  that  for  the  important  paths.  .  . 

(w(z)-w(z'))2  +  (z-z ' )  2  (z-z ' )  2  (76) 

[in  our  notation  for  paths]."  It  is  not  hard  to  see  that  (76)  and  (73)  are 
essentially  equivalent.  However,  here  the  approximation  is  explicitly  made  only 
for  "important  paths."  It  is  clear  from  Dashen 's  quoted  explanation  that  he  has 
in  mind  the  physically  important  rays  in  situations  where  the  small-angle 
approximation  leading  to  the  Parabolic  Wave  Equation  is  valid.  But  there  seems 
to  be  no  clear  mathematical  basis  for  neglecting  paths  that  wander  away  from  the 
z-axis  in  Dashen's  Feynman  integrals:  the  "measure"  in  path  space  corresponding 
to  Feynman  integrals,  if  we  could  construct  one,  would  be  infinite-dimensional 
Lebesque  measure,  and  therefore  translation  invariant.  Much  the  same  criticism 
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would  seem  to  apply  to  the  Brownian  paths  in  our  Wiener  integral  formula  (64) . 

In  order  to  make  the  Markov  approximation  for  "important  paths"  more  nearly 
exact,  we  should  use  a  measure  in  path  space  such  that  paths  which  wander  far 
from  the  z-axis  are  truly  negligible,  i.e.  constitute  a  set  of  negligibly  small 
measure.  An  ideal  candidate  for  such  a  measure  would  seem  to  be  the  measure 
corresponding  to  the  random  ray  process  of  section  III.  Assuming  that  the 
parabolic  approximation  is  self-consistent,  the  quantity  | ip (x , z)  |  2 ,  which  is 
proportional  to  the  acoustic  intensity,  should  be  concentrated  near  the  z-axis. 
But  |ty(x,z) |  =  p(x,z) ,  the  probability  density  of  the  position  x(z)  of  the 

random  ray  at  range  z.  Thus,  in  the  sequel  we  shall  use  the  Cameron-Martin- 
Girsanov  formula  ([12];  see  also  [30,  chapter  7],  where  this  name  is  not  used) 
to  change  variables  in  the  path  integral  of  equation  (64)  from  the  Brownian  path 
w(s)  to  the  random  ray  x  (s)  . 
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